Load packages

# rm(list = ls())
library(tictoc)  # to monitor time

library(raster) # read raster files
library(rgdal) # Use GDAL functions
library(exactextractr) # Zonal statistics from raster
library(sf) # Working with shapefile layers

library(tidyverse) # data wrangling
library(plyr) # data wrangling
library(dplyr) # data wrangling

library(ggpmisc) # To place geom_text within plot
library(plotly) # Interactive visualization
library(patchwork) # This package for arranging multiple ggplots

print("Packages loaded successfully!")
## [1] "Packages loaded successfully!"

Functions specific to this RMD

#########################################
# Creating function for exponential fit
myexpfit <- function(response_var, pred_var, mydata){
  # print(paste0("Model = ", "log(",response_var, ")~", pred_var))
  model.linear <- lm(paste0("log(",response_var, ")~", pred_var), data = mydata)
   a1 <- exp(model.linear$coefficients[1])
   b1 <- (model.linear$coefficients[2])
   R2 <- summary(model.linear)$r.squared

   myresults <- c(a1, b1, R2)
   return(myresults)
}
#########################################

# My GGPLOT for yield map
my_yield_plot <- function(yield_df){
  fieldtitle <- levels(yield_df$Fieldname)
  # Legend labeling - Splitting the yield range into equal classes
  minyld <- round(min(yield_df$Yield_Mg.ha),0) # minimum yield 
  maxyld <- round(max(yield_df$Yield_Mg.ha),0) # maximum yield
  Nclasses <- 4      # Number of yield classes to split from the whole yield range
  
  mybreaks <- round(seq(minyld, maxyld, len = (Nclasses+1)),0)
  YE_colors <-  c("blue","cyan","lightgreen","yellow","orange", "red")

  yld_plot <- ggplot(yield_df, aes(x = Longitude, y = Latitude)) +
  geom_point(aes(color = Yield_Mg.ha), size = 0.7) +    
      theme_bw() +
      # facet_grid(Field~Year, scales = "free") +
      theme(aspect.ratio = 1.0) +
      # ggtitle("Actual yield (Mg/ha)")+
      facet_wrap(~Fieldname, scales = "free") +
      # theme(legend.position = c(0.85, 0.27))+
      scale_colour_gradientn(colours = YE_colors,  
                             name = "Yield (Mg/ha)", 
                             limits=c(minyld,maxyld),
                             breaks = mybreaks,
                             labels = round(mybreaks, 0)) +
    guides(colour = guide_colourbar(barheight = unit(1.8, "in")))+
    theme(panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(), 
            axis.title.x = element_blank(),
            axis.title.y = element_blank(),
            axis.text.x = element_blank(),
            axis.text.y = element_blank(),
            axis.ticks = element_blank(), 
            plot.title = element_text(size=14, face = "bold")) +
      theme(legend.title = element_text(face = "bold", size=8, margin = margin(0,0,7,0)),
            legend.text = element_text(vjust = 0.5),
            strip.text = element_text(face = "bold", size = 10))
  yld_plot
  
}

#########################################
my_pred_plot <- function(sub_dat){
  fieldtitle <- levels(sub_dat$Fieldname)
    # Legend labeling - Splitting the yield range into equal classes
  minyld <- round(min(sub_dat$Pred_Yield),0) # minimum yield 
  maxyld <- round(max(sub_dat$Pred_Yield),0) # maximum yield
  Nclasses <- 4      # Number of yield classes to split from the whole yield range
    
  mybreaks <- round(seq(minyld, maxyld, len = (Nclasses+1)),0)
  YE_colors <-  c("blue","cyan","lightgreen","yellow","orange", "red")

  pred_map <- ggplot(sub_dat, aes(x = Longitude, y = Latitude)) +
        geom_point(aes(color = Pred_Yield), size = 0.7) +    
        theme_bw() +
        # facet_grid(Field~Year, scales = "free") +
        theme(aspect.ratio = 1.0) +
        # ggtitle("Predicted yield (Mg/ha)")+
        facet_wrap(~Fieldname, scales = "free") +
        # theme(legend.position = c(0.85, 0.27))+
        scale_colour_gradientn(colours = YE_colors,  
                               name = "Yield (Mg/ha)", 
                               limits=c(minyld,maxyld),
                               breaks = mybreaks,
                               labels = round(mybreaks, 0)) +
        guides(colour = guide_colourbar(barheight = unit(1.8, "in")))+
        theme(panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(), 
              axis.title.x = element_blank(),
              axis.title.y = element_blank(),
              axis.text.x = element_blank(),
              axis.text.y = element_blank(),
              axis.ticks = element_blank(), 
              plot.title = element_text(size=14, face = "bold")) +
        theme(legend.title = element_text(face = "bold", size=8, margin = margin(0,0,7,0)),
              legend.text = element_text(vjust = 0.5),
              strip.text = element_text(face = "bold", size = 10))
  pred_map
}

#########################################
## Actual and predicted yield map
make_yield_maps <- function(sub_datf){

  sel_dat4 <- sub_datf %>% 
    select(Fieldname, Latitude, Longitude, Yield_Mg.ha, Pred_Yield) %>% 
    gather("yield_type", "Value", 4:5) %>% 
    droplevels()
  
  sel_dat4$yield_type <- recode_factor(factor(sel_dat4$yield_type), Pred_Yield = "Predicted yield", Yield_Mg.ha = "Actual yield")
  
  sel_dat4$yield_type <- factor(sel_dat4$yield_type, levels = c("Actual yield", "Predicted yield"))
  
  minyld <- round(min(sel_dat4$Value),0) # minimum yield 
  maxyld <- round(max(sel_dat4$Value),0) # maximum yield
  Nclasses <- 4      # Number of yield classes to split from the whole yield range
    
  mybreaks <- round(seq(minyld, maxyld, len = (Nclasses+1)),0)
  YE_colors <-  c("blue","cyan","lightgreen","yellow","orange", "red")
  
  yld_plot <- ggplot(sel_dat4, aes(x = Longitude, y = Latitude)) +
    geom_point(aes(color = Value), size = 0.7) +    
        theme_bw() +
        # facet_grid(Field~Year, scales = "free") +
        theme(aspect.ratio = 1.0) +
        # ggtitle("Actual yield (Mg/ha)")+
        facet_wrap(~yield_type, scales = "free") +
        # theme(legend.position = c(0.85, 0.27))+
        scale_colour_gradientn(colours = YE_colors,  
                               name = "Yield (Mg/ha)", 
                               limits=c(minyld,maxyld),
                               breaks = mybreaks,
                               labels = round(mybreaks, 0)) +
      guides(colour = guide_colourbar(barheight = unit(1.8, "in")))+
      theme(panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(), 
              axis.title.x = element_blank(),
              axis.title.y = element_blank(),
              axis.text.x = element_blank(),
              axis.text.y = element_blank(),
              axis.ticks = element_blank(), 
              plot.title = element_text(size=14, face = "bold")) +
        theme(legend.title = element_text(face = "bold", size=8, margin = margin(0,0,7,0)),
              legend.text = element_text(vjust = 0.5),
              strip.text = element_text(face = "bold", size = 10))
  
  yld_plot
}

Load dataset

load("Quantix_deep_learning.RData", verbose = TRUE)
## Loading objects:
##   dat
str(dat)
## 'data.frame':    1270534 obs. of  46 variables:
##  $ Fieldname            : Factor w/ 8 levels "DMD_GH1","PSF_111",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ FlightDate           : Factor w/ 22 levels "2019/06/07","2019/06/09",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ DOY                  : int  158 158 158 158 158 158 158 158 158 158 ...
##  $ Week                 : Factor w/ 14 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Latitude             : num  43 43 43 43 43 ...
##  $ Longitude            : num  -76.6 -76.6 -76.6 -76.6 -76.6 ...
##  $ Yield                : num  182 188 191 194 197 ...
##  $ Yield_Mg.ha          : num  11.4 11.8 12 12.2 12.4 ...
##  $ Type                 : Factor w/ 2 levels "grain","silage": 1 1 1 1 1 1 1 1 1 1 ...
##  $ N                    : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Strip                : Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ rgn_red              : num  18 18 18 18 18 18 18 18 18 18 ...
##  $ rgn_green            : num  16 15 15 15 15 15 15 15 15 16 ...
##  $ rgn_nir              : num  36 36 35 35 35 35 35 35 35 35 ...
##  $ rgb_red              : num  115 115 116 117 118 118 119 118 118 119 ...
##  $ rgb_green            : num  89 89 90 90 91 91 92 91 91 92 ...
##  $ rgb_blue             : num  67 67 68 68 69 69 70 69 69 70 ...
##  $ NDVI                 : num  0.333 0.333 0.321 0.321 0.321 ...
##  $ GNDVI                : num  0.385 0.412 0.4 0.4 0.4 ...
##  $ EVI2                 : num  0.597 0.616 0.59 0.59 0.59 ...
##  $ SR                   : num  2 2 1.94 1.94 1.94 ...
##  $ EXG                  : num  -4 -4 -4 -5 -5 -5 -5 -5 -5 -5 ...
##  $ TGI                  : num  3.28 3.28 3.28 2.89 2.89 2.89 2.89 2.89 2.89 2.89 ...
##  $ Mgmt_zone            : Factor w/ 4 levels "1","2","3","4": 3 3 3 3 3 3 2 2 2 2 ...
##  $ Climod_GDD           : int  10 10 10 10 10 10 10 10 10 10 ...
##  $ Climod_Temp_C        : num  15.6 15.6 15.6 15.6 15.6 ...
##  $ Climod_Precip_mm     : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ ACIS_GDD             : int  10 10 10 10 10 10 10 10 10 10 ...
##  $ ACIS_Temp_C          : num  15.6 15.6 15.6 15.6 15.6 ...
##  $ ACIS_Precip_mm       : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ NASA_PAR             : num  159 159 159 159 159 ...
##  $ NASA_Temp_C          : num  17.3 17.3 17.3 17.3 17.3 ...
##  $ NASA_SH              : num  7.2 7.2 7.2 7.2 7.2 7.2 7.2 7.2 7.2 7.2 ...
##  $ NASA_RH              : num  59.7 59.7 59.7 59.7 59.7 ...
##  $ NASA_Precip_mm       : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ NASA_WS              : num  0.56 0.56 0.56 0.56 0.56 0.56 0.56 0.56 0.56 0.56 ...
##  $ NASA_SurfWet         : num  0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 ...
##  $ NASA_RootWet         : num  0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 ...
##  $ NASA_ProfWet         : num  0.66 0.66 0.66 0.66 0.66 0.66 0.66 0.66 0.66 0.66 ...
##  $ DEM                  : int  125 125 125 125 125 125 125 125 125 125 ...
##  $ DEM_type_Flat_Area   : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ DEM_type_Lower_Slope : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 2 2 2 ...
##  $ DEM_type_Middle_Slope: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ DEM_type_Ridge       : Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 1 1 1 ...
##  $ DEM_type_Upper_Slope : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ DEM_type_Valley      : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
r2df_all <- data.frame()

Data preprocessing

Check for missing values

# Check if any columns have missing values 
sapply(dat, anyNA)
##             Fieldname            FlightDate                   DOY 
##                 FALSE                 FALSE                 FALSE 
##                  Week              Latitude             Longitude 
##                 FALSE                 FALSE                 FALSE 
##                 Yield           Yield_Mg.ha                  Type 
##                 FALSE                 FALSE                 FALSE 
##                     N                 Strip               rgn_red 
##                  TRUE                  TRUE                 FALSE 
##             rgn_green               rgn_nir               rgb_red 
##                 FALSE                 FALSE                 FALSE 
##             rgb_green              rgb_blue                  NDVI 
##                 FALSE                 FALSE                 FALSE 
##                 GNDVI                  EVI2                    SR 
##                 FALSE                 FALSE                 FALSE 
##                   EXG                   TGI             Mgmt_zone 
##                 FALSE                 FALSE                  TRUE 
##            Climod_GDD         Climod_Temp_C      Climod_Precip_mm 
##                 FALSE                 FALSE                 FALSE 
##              ACIS_GDD           ACIS_Temp_C        ACIS_Precip_mm 
##                 FALSE                 FALSE                 FALSE 
##              NASA_PAR           NASA_Temp_C               NASA_SH 
##                 FALSE                 FALSE                 FALSE 
##               NASA_RH        NASA_Precip_mm               NASA_WS 
##                 FALSE                 FALSE                 FALSE 
##          NASA_SurfWet          NASA_RootWet          NASA_ProfWet 
##                 FALSE                 FALSE                 FALSE 
##                   DEM    DEM_type_Flat_Area  DEM_type_Lower_Slope 
##                 FALSE                 FALSE                 FALSE 
## DEM_type_Middle_Slope        DEM_type_Ridge  DEM_type_Upper_Slope 
##                 FALSE                 FALSE                 FALSE 
##       DEM_type_Valley 
##                 FALSE
# Names of columns that has missing values
names(dat)[sapply(dat, anyNA)]
## [1] "N"         "Strip"     "Mgmt_zone"

Plot all yield data

g_dat <- dat %>%
  dplyr::filter(Type == "grain") %>%
  droplevels()

my_yield_plot(g_dat)

s_dat <- dat %>%
  dplyr::filter(Type == "silage") %>%
  droplevels()

my_yield_plot(s_dat)

Exploratory data analysis

Distribution plots

# Violin plot

vplt <- ggplot(dat, aes(x = Fieldname, y = Yield_Mg.ha, fill = Fieldname)) +
  geom_violin(alpha = 0.4) +
  theme_bw() +
  facet_wrap(~Type, nrow = 2, scales = "free_y")

vplt

dplt <- ggplot(dat, aes(x = Yield_Mg.ha, y = ..scaled.., fill = Fieldname)) +
  geom_density(alpha = 0.4, size = 1) + 
  theme_bw() +
  facet_wrap(~Type, nrow = 2, scales = "free")

dplt

GPS points vs yield

# "DMD_GH1" "PSF_111" "PSF_12"  "SLS_ABH" "SLS_NS"  "SSF_121" "SSF_202" "SSF_66"
# selfield <- "DMD_GH1"

sub_datf <- dat %>% 
  filter(Week == 1 & Type == "grain") %>% 
  droplevels()

# Legend labeling - Splitting the yield range into equal classes
minyld <- round(min(sub_datf$Yield_Mg.ha),0) # minimum yield 
maxyld <- round(max(sub_datf$Yield_Mg.ha),0) # maximum yield
Nclasses <- 4      # Number of yield classes to split from the whole yield range

mybreaks <- round(seq(minyld, maxyld, len = (Nclasses+1)),0)
YE_colors <-  c("blue","cyan","lightgreen","yellow","orange", "red")

#=====================
lat_plot <- ggplot(sub_datf, aes(x = Latitude, y = Yield_Mg.ha)) +
  geom_point(aes(color = Yield_Mg.ha), size = 0.7) +    
      theme_bw() +
      # facet_grid(Field~Year, scales = "free") +
      theme(aspect.ratio = 1.0) +
      ggtitle("Latitude vs Yield")+
      facet_wrap(~Fieldname, scales = "free") +
      # theme(legend.position = c(0.85, 0.27))+
      scale_colour_gradientn(colours = YE_colors,  
                             name = "Yield (Mg/ha)", 
                             limits=c(minyld,maxyld),
                             breaks = mybreaks,
                             labels = round(mybreaks, 0)) +
    guides(colour = guide_colourbar(barheight = unit(1.8, "in")))+
    theme(panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(), 
            axis.text.x = element_blank(),
            axis.text.y = element_blank(),
            axis.ticks = element_blank(), 
            plot.title = element_text(size=14, face = "bold")) +
      theme(legend.title = element_text(face = "bold", size=8, margin = margin(0,0,7,0)),
            legend.text = element_text(vjust = 0.5),
            strip.text = element_text(face = "bold", size = 10))
lat_plot

#=====================
long_plot <- ggplot(sub_datf, aes(x = Longitude, y = Yield_Mg.ha)) +
  geom_point(aes(color = Yield_Mg.ha), size = 0.7) +    
      theme_bw() +
      # facet_grid(Field~Year, scales = "free") +
      theme(aspect.ratio = 1.0) +
      ggtitle("Longitude vs Yield")+
      facet_wrap(~Fieldname, scales = "free") +
      # theme(legend.position = c(0.85, 0.27))+
      scale_colour_gradientn(colours = YE_colors,  
                             name = "Yield (Mg/ha)", 
                             limits=c(minyld,maxyld),
                             breaks = mybreaks,
                             labels = round(mybreaks, 0)) +
    guides(colour = guide_colourbar(barheight = unit(1.8, "in")))+
    theme(panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(), 
            axis.text.x = element_blank(),
            axis.text.y = element_blank(),
            axis.ticks = element_blank(), 
            plot.title = element_text(size=14, face = "bold")) +
      theme(legend.title = element_text(face = "bold", size=8, margin = margin(0,0,7,0)),
            legend.text = element_text(vjust = 0.5),
            strip.text = element_text(face = "bold", size = 10))
long_plot

Approach 1

1.a. Use only strip data for exponential model

# "DMD_GH1" "PSF_111" "PSF_12"  "SLS_ABH" "SLS_NS"  "SSF_121" "SSF_202" "SSF_66"
# "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14"
# Best grain field - DMD_GH1
# Best silage filed - SSF_121, PSF_12

vi <- c("NDVI", "GNDVI", "EVI2", "SR", "EXG", "TGI")

selfield <- "DMD_GH1"
# seltype <- "silage"

sub_dat <- dat %>% 
  filter(!is.na(N) & Fieldname == selfield) %>% 
  droplevels()

my_yield_plot(sub_dat)

# "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14"
lvl <- levels(sub_dat$Week)
sub_dat$Strip <- as.numeric(sub_dat$Strip)

#==============
# Scatter plot between NDVI and Yield
sct_plt <- ggplot(sub_dat, aes(x = NDVI, y = Yield_Mg.ha)) + 
  geom_point(aes(color = Week)) +
  theme_bw() + 
  theme(aspect.ratio = 1) +
  ggtitle("NDVI vs Yield") +
  xlab("NDVI") + 
  ylab("Yield (Mg/ha)") + 
  xlim(c(0.0, 1.0))

sct_plt

#==============

r2df <- data.frame()
for (i in 1:length(lvl)){
  # print(i)
  wsub_dat <- sub_dat %>% 
    filter(Week == lvl[i]) %>%
    droplevels()
  
  for (j in 1:length(vi)){
      
    dep_var <- "Yield_Mg.ha"
    ind_var <- vi[j]
    # cat(dep_var, "---", ind_var, "\n")
    
    mod <- myexpfit("Yield_Mg.ha", ind_var, wsub_dat)
    
    ## RMSE calculation ======
    a_fit <- as.numeric(mod[1])
    b_fit <- as.numeric(mod[2])
    rsq <- as.numeric(mod[3])
    
    ssub_dat <- dat %>% 
      filter(Week == lvl[i] & Fieldname == selfield) %>% 
      droplevels() 
      
    ssub_dat$Pred_Yield <-  a_fit*exp(b_fit*ssub_dat[[ind_var]])
  
    ssub_dat$error <- (ssub_dat$Pred_Yield - ssub_dat$Yield_Mg.ha)
    sqe <- (ssub_dat$error^2)
    mse <- mean(sqe)
    rmse_week <- sqrt(mse)
    # cat(lvl[i], ") a_fit = ", a_fit, "; b_fit = ", b_fit, ";\t R2 = ", rsq,";\t RMSE = ", rmse_week, "\n")
    # cat("RMSE for", selfield, "on week", lvl[i], "is", rmse, "Mg/ha")
    ## ========================
    
    r2 <- c(as.numeric(c(lvl[i], mod[1], mod[2], mod[3], rmse_week)), ind_var)
    r2df <- rbind(r2df, r2)
  }

  # print(paste0(i, " = ", mod[3]))
}

colnames(r2df) <- c("Weeks", "a", "b", "R2", "RMSE", "VI")

num_cols <- c("Weeks", "a", "b", "R2", "RMSE")
r2df[num_cols] <- lapply(r2df[num_cols], as.numeric)
str(r2df)
## 'data.frame':    78 obs. of  6 variables:
##  $ Weeks: num  1 1 1 1 1 1 2 2 2 2 ...
##  $ a    : num  11.13 10.62 11.21 9.89 14.35 ...
##  $ b    : num  0.87957 0.84782 0.47224 0.20573 -0.00694 ...
##  $ R2   : num  0.0329 0.0314 0.0371 0.0349 0.0249 ...
##  $ RMSE : num  2.29 2.21 2.3 3.85 2.02 ...
##  $ VI   : chr  "NDVI" "GNDVI" "EVI2" "SR" ...
r2df$Fieldname <- rep(selfield, length(lvl))
r2df$Approach <- rep("1a", length(lvl))
r2df_all <- rbind(r2df_all, r2df)

1.b. Use lower and upper 10% of strip data

# "DMD_GH1" "PSF_111" "PSF_12"  "SLS_ABH" "SLS_NS"  "SSF_121" "SSF_202" "SSF_66"
selfield <- "SSF_66"

sub_dat <- dat %>% 
  filter(!is.na(N) & Fieldname == selfield) %>%
  droplevels()
my_yield_plot(sub_dat)

# Evaluate qunatile cutoffs
pert <- 0.10

cutoffs <- quantile(sub_dat$Yield_Mg.ha, c(pert, (1-pert)))
cat("Cutoff for field", selfield, " =", cutoffs)
## Cutoff for field SSF_66  = 59.6752 81.1395
# Subset points within cutoff
qsub_dat <- sub_dat %>% 
  filter(Yield_Mg.ha <= cutoffs[[1]] | Yield_Mg.ha >= cutoffs[[2]]) %>% 
  droplevels()

my_yield_plot(qsub_dat)

#==============
# Scatter plot between NDVI and Yield
sct_plt <- ggplot(qsub_dat, aes(x = NDVI, y = Yield_Mg.ha)) + 
  geom_point(aes(color = Week)) +
  theme_bw() + 
  theme(aspect.ratio = 1) +
  xlab("NDVI") + 
  ylab("Yield (Mg/ha)") + 
  xlim(c(0.0, 1.0))

sct_plt

#==============

lvl <- levels(qsub_dat$Week)

r2df <- data.frame()
for (i in 1:length(lvl)){
  # print(i)
  wsub_dat <- qsub_dat %>% 
    filter(Week == lvl[i]) %>%
    droplevels()

    for (j in 1:length(vi)){
      
      dep_var <- "Yield_Mg.ha"
      ind_var <- vi[j]
      # cat(dep_var, "---", ind_var, "\n")
      
      mod <- myexpfit("Yield_Mg.ha", ind_var, wsub_dat)
     
      ## RMSE calculation ======
      a_fit <- as.numeric(mod[1])
      b_fit <- as.numeric(mod[2])
      rsq <- as.numeric(mod[3])
      
      ssub_dat <- dat %>% 
        filter(Week == lvl[i] & Fieldname == selfield) %>% 
        droplevels() 
        
      ssub_dat$Pred_Yield <-  a_fit*exp(b_fit*ssub_dat[[ind_var]])
    
      ssub_dat$error <- (ssub_dat$Pred_Yield - ssub_dat$Yield_Mg.ha)
      sqe <- (ssub_dat$error^2)
      mse <- mean(sqe)
      rmse_week <- sqrt(mse)
      # cat(lvl[i], ") a_fit = ", a_fit, "; b_fit = ", b_fit, ";\t R2 = ", rsq,";\t RMSE = ", rmse_week, "\n")
      # cat("RMSE for", selfield, "on week", lvl[i], "is", rmse, "Mg/ha")
      ## ========================
      
      r2 <- c(as.numeric(c(lvl[i], mod[1], mod[2], mod[3], rmse_week)), ind_var)
      r2df <- rbind(r2df, r2)
    }
  # print(paste0(i, " = ", mod[3]))
}

colnames(r2df) <- c("Weeks", "a", "b", "R2", "RMSE", "VI")

num_cols <- c("Weeks", "a", "b", "R2", "RMSE")
r2df[num_cols] <- lapply(r2df[num_cols], as.numeric)
str(r2df)
## 'data.frame':    78 obs. of  6 variables:
##  $ Weeks: num  1 1 1 1 1 1 2 2 2 2 ...
##  $ a    : num  82.6 101.3 83.4 80.3 67.1 ...
##  $ b    : num  -0.64658 -1.11741 -0.38113 -0.08987 0.00485 ...
##  $ R2   : num  0.0106 0.0228 0.01513 0.01428 0.00636 ...
##  $ RMSE : num  8.97 9 8.97 8.97 9.16 ...
##  $ VI   : chr  "NDVI" "GNDVI" "EVI2" "SR" ...
r2df$Fieldname <- rep(selfield, length(lvl))
r2df$Approach <- rep("1b", length(lvl))
r2df_all <- rbind(r2df_all, r2df)

Approach 2

2.a. Use whole field data without strip –> is.na(N)

# "DMD_GH1" "PSF_111" "PSF_12"  "SLS_ABH" "SLS_NS"  "SSF_121" "SSF_202" "SSF_66"
selfield <- "SSF_66"

sub_dat <- dat %>% 
  filter(is.na(N) & Fieldname == selfield) %>% 
  droplevels()

my_yield_plot(sub_dat)

#==============
# Scatter plot between NDVI and Yield
sct_plt <- ggplot(sub_dat, aes(x = NDVI, y = Yield_Mg.ha)) + 
  geom_point(aes(color = Week)) +
  theme_bw() + 
  theme(aspect.ratio = 1) +
  xlab("NDVI") + 
  ylab("Yield (Mg/ha)") + 
  xlim(c(0.0, 1.0))

sct_plt

#==============

lvl <- levels(sub_dat$Week)

r2df <- data.frame()
for (i in 1:length(lvl)){
  # print(i)
  wsub_dat <- sub_dat %>% 
    filter(Week == lvl[i]) %>%
    droplevels()
  
    for (j in 1:length(vi)){
      
      dep_var <- "Yield_Mg.ha"
      ind_var <- vi[j]
      # cat(dep_var, "---", ind_var, "\n")
      
      mod <- myexpfit("Yield_Mg.ha", ind_var, wsub_dat)
    
      ## RMSE calculation ======
      a_fit <- as.numeric(mod[1])
      b_fit <- as.numeric(mod[2])
      rsq <- as.numeric(mod[3])
      
      ssub_dat <- dat %>% 
        filter(Week == lvl[i] & Fieldname == selfield) %>% 
        droplevels() 
        
      ssub_dat$Pred_Yield <-  a_fit*exp(b_fit*ssub_dat[[ind_var]])
    
      ssub_dat$error <- (ssub_dat$Pred_Yield - ssub_dat$Yield_Mg.ha)
      sqe <- (ssub_dat$error^2)
      mse <- mean(sqe)
      rmse_week <- sqrt(mse)
      # cat(lvl[i], ") a_fit = ", a_fit, "; b_fit = ", b_fit, ";\t R2 = ", rsq,";\t RMSE = ", rmse_week, "\n")
      # cat("RMSE for", selfield, "on week", lvl[i], "is", rmse, "Mg/ha")
      ## ========================
      
      r2 <- c(as.numeric(c(lvl[i], mod[1], mod[2], mod[3], rmse_week)), ind_var)
      r2df <- rbind(r2df, r2)
  }
  # print(paste0(i, " = ", mod[3]))
}

colnames(r2df) <- c("Weeks", "a", "b", "R2", "RMSE", "VI")

num_cols <- c("Weeks", "a", "b", "R2", "RMSE")
r2df[num_cols] <- lapply(r2df[num_cols], as.numeric)
str(r2df)
## 'data.frame':    78 obs. of  6 variables:
##  $ Weeks: num  1 1 1 1 1 1 2 2 2 2 ...
##  $ a    : num  94.2 101 90.9 85.6 68.5 ...
##  $ b    : num  -1.06597 -1.09442 -0.5315 -0.12254 -0.00442 ...
##  $ R2   : num  0.0535 0.0397 0.0523 0.0393 0.0132 ...
##  $ RMSE : num  8.96 8.96 8.95 8.94 9.09 ...
##  $ VI   : chr  "NDVI" "GNDVI" "EVI2" "SR" ...
r2df$Fieldname <- rep(selfield, length(lvl))
r2df$Approach <- rep("2a", length(lvl))
r2df_all <- rbind(r2df_all, r2df)

2.b. Use lower and upper 10% of whole field data

# "DMD_GH1" "PSF_111" "PSF_12"  "SLS_ABH" "SLS_NS"  "SSF_121" "SSF_202" "SSF_66"
selfield <- "SSF_66"

sub_dat <- dat %>% 
  filter(is.na(N) & Fieldname == selfield) %>% 
  droplevels()

my_yield_plot(sub_dat)

# Evaluate quantile cutoffs
pert <- 0.10
cutoffs <- quantile(sub_dat$Yield_Mg.ha, c(pert, (1-pert)))
cat("Cutoff for field", selfield, " =", cutoffs)
## Cutoff for field SSF_66  = 56.8841 81.4359
# Subset points within cutoff
qsub_dat <- sub_dat %>% 
  filter(Yield_Mg.ha <= cutoffs[[1]] | Yield_Mg.ha >= cutoffs[[2]]) %>% 
  droplevels()

my_yield_plot(qsub_dat)

#==============
# Scatter plot between NDVI and Yield
sct_plt <- ggplot(qsub_dat, aes(x = NDVI, y = Yield_Mg.ha)) + 
  geom_point(aes(color = Week)) +
  theme_bw() + 
  theme(aspect.ratio = 1) +
  xlab("NDVI") + 
  ylab("Yield (Mg/ha)") + 
  xlim(c(0.0, 1.0))

sct_plt

#==============

lvl <- levels(qsub_dat$Week)

r2df <- data.frame()
for (i in 1:length(lvl)){
  # print(i)
  wsub_dat <- qsub_dat %>% 
    filter(Week == lvl[i]) %>%
    droplevels()
  
    for (j in 1:length(vi)){
      
      dep_var <- "Yield_Mg.ha"
      ind_var <- vi[j]
      # cat(dep_var, "---", ind_var, "\n")
      
      mod <- myexpfit("Yield_Mg.ha", ind_var, wsub_dat)
      
      ## RMSE calculation ======
      a_fit <- as.numeric(mod[1])
      b_fit <- as.numeric(mod[2])
      rsq <- as.numeric(mod[3])
      
      ssub_dat <- dat %>% 
        filter(Week == lvl[i] & Fieldname == selfield) %>% 
        droplevels() 
        
      ssub_dat$Pred_Yield <-  a_fit*exp(b_fit*ssub_dat[[ind_var]])
    
      ssub_dat$error <- (ssub_dat$Pred_Yield - ssub_dat$Yield_Mg.ha)
      sqe <- (ssub_dat$error^2)
      mse <- mean(sqe)
      rmse_week <- sqrt(mse)
      cat(lvl[i], ") a_fit = ", a_fit, "; b_fit = ", b_fit, ";\t R2 = ", rsq,";\t RMSE = ", rmse_week, "\n")
      # cat("RMSE for", selfield, "on week", lvl[i], "is", rmse, "Mg/ha")
      ## ========================
  
      r2 <- c(as.numeric(c(lvl[i], mod[1], mod[2], mod[3], rmse_week)), ind_var)  
      r2df <- rbind(r2df, r2)
  }
  # print(paste0(i, " = ", mod[3]))
}
## 1 ) a_fit =  114.5396 ; b_fit =  -1.775352 ;  R2 =  0.1047697 ;   RMSE =  9.3737 
## 1 ) a_fit =  127.595 ; b_fit =  -1.814529 ;   R2 =  0.07242787 ;  RMSE =  9.418363 
## 1 ) a_fit =  104.0244 ; b_fit =  -0.8232189 ;     R2 =  0.0940588 ;   RMSE =  9.351495 
## 1 ) a_fit =  85.56061 ; b_fit =  -0.138763 ;  R2 =  0.05338813 ;  RMSE =  9.445993 
## 1 ) a_fit =  69.44581 ; b_fit =  -0.01730106 ;    R2 =  0.08496018 ;  RMSE =  9.968001 
## 1 ) a_fit =  80.12686 ; b_fit =  -0.02549843 ;    R2 =  0.0559519 ;   RMSE =  9.91363 
## 2 ) a_fit =  97.68063 ; b_fit =  -0.8558484 ;     R2 =  0.08915737 ;  RMSE =  9.710312 
## 2 ) a_fit =  100.4366 ; b_fit =  -0.8612361 ;     R2 =  0.05955035 ;  RMSE =  9.634497 
## 2 ) a_fit =  90.54516 ; b_fit =  -0.3716412 ;     R2 =  0.08503231 ;  RMSE =  9.616456 
## 2 ) a_fit =  75.43675 ; b_fit =  -0.04808592 ;    R2 =  0.08371296 ;  RMSE =  9.359314 
## 2 ) a_fit =  72.97827 ; b_fit =  -0.01113459 ;    R2 =  0.07962387 ;  RMSE =  10.06643 
## 2 ) a_fit =  77.92035 ; b_fit =  -0.01842857 ;    R2 =  0.05551291 ;  RMSE =  9.937669 
## 3 ) a_fit =  130.2284 ; b_fit =  -1.772073 ;  R2 =  0.08667496 ;  RMSE =  9.545287 
## 3 ) a_fit =  151.7014 ; b_fit =  -1.912627 ;  R2 =  0.07356775 ;  RMSE =  9.665523 
## 3 ) a_fit =  120.4115 ; b_fit =  -0.8546063 ;     R2 =  0.08906408 ;  RMSE =  9.555316 
## 3 ) a_fit =  93.15247 ; b_fit =  -0.1546254 ;     R2 =  0.05405708 ;  RMSE =  9.430931 
## 3 ) a_fit =  77.07453 ; b_fit =  -0.02416504 ;    R2 =  0.07798439 ;  RMSE =  9.901467 
## 3 ) a_fit =  87.24738 ; b_fit =  -0.03597171 ;    R2 =  0.05820835 ;  RMSE =  10.1733 
## 4 ) a_fit =  39.96259 ; b_fit =  0.7992891 ;  R2 =  0.049653 ;    RMSE =  9.760227 
## 4 ) a_fit =  35.31187 ; b_fit =  0.9765932 ;  R2 =  0.05753524 ;  RMSE =  9.731614 
## 4 ) a_fit =  44.29221 ; b_fit =  0.3180522 ;  R2 =  0.04569029 ;  RMSE =  9.737102 
## 4 ) a_fit =  55.23613 ; b_fit =  0.03822832 ;     R2 =  0.0205722 ;   RMSE =  9.777046 
## 4 ) a_fit =  80.0245 ; b_fit =  -0.01120671 ;     R2 =  0.02603877 ;  RMSE =  9.936435 
## 4 ) a_fit =  95.7974 ; b_fit =  -0.03033375 ;     R2 =  0.05545126 ;  RMSE =  9.916124 
## 5 ) a_fit =  21.26627 ; b_fit =  1.377981 ;   R2 =  0.1782023 ;   RMSE =  9.033781 
## 5 ) a_fit =  20.33525 ; b_fit =  1.457281 ;   R2 =  0.1927672 ;   RMSE =  9.040896 
## 5 ) a_fit =  26.53846 ; b_fit =  0.5126518 ;  R2 =  0.2002983 ;   RMSE =  9.111857 
## 5 ) a_fit =  42.33944 ; b_fit =  0.04008052 ;     R2 =  0.1821739 ;   RMSE =  9.767657 
## 5 ) a_fit =  50.62741 ; b_fit =  0.007041083 ;    R2 =  0.05393396 ;  RMSE =  9.303907 
## 5 ) a_fit =  52.54574 ; b_fit =  0.01057815 ;     R2 =  0.04159836 ;  RMSE =  9.362912 
## 6 ) a_fit =  17.90351 ; b_fit =  1.559624 ;   R2 =  0.153909 ;    RMSE =  9.08576 
## 6 ) a_fit =  15.60691 ; b_fit =  1.747906 ;   R2 =  0.163055 ;    RMSE =  9.040844 
## 6 ) a_fit =  23.14956 ; b_fit =  0.5700638 ;  R2 =  0.1629603 ;   RMSE =  9.067256 
## 6 ) a_fit =  40.48376 ; b_fit =  0.04118642 ;     R2 =  0.1547518 ;   RMSE =  9.504802 
## 6 ) a_fit =  75.83658 ; b_fit =  -0.004390261 ;   R2 =  0.009154759 ;     RMSE =  9.77616 
## 6 ) a_fit =  84.54597 ; b_fit =  -0.01349674 ;    R2 =  0.02717307 ;  RMSE =  9.746971 
## 7 ) a_fit =  12.22936 ; b_fit =  1.892968 ;   R2 =  0.1274779 ;   RMSE =  9.175553 
## 7 ) a_fit =  10.66996 ; b_fit =  2.097191 ;   R2 =  0.1527872 ;   RMSE =  9.029418 
## 7 ) a_fit =  18.26674 ; b_fit =  0.6438008 ;  R2 =  0.1446185 ;   RMSE =  9.052002 
## 7 ) a_fit =  43.57484 ; b_fit =  0.0227518 ;  R2 =  0.1088754 ;   RMSE =  9.814816 
## 7 ) a_fit =  84.75155 ; b_fit =  -0.007686897 ;   R2 =  0.03953933 ;  RMSE =  10.2967 
## 7 ) a_fit =  85.4768 ; b_fit =  -0.01445088 ;     R2 =  0.04640809 ;  RMSE =  10.3253 
## 8 ) a_fit =  7.357011 ; b_fit =  2.434166 ;   R2 =  0.1423155 ;   RMSE =  8.993199 
## 8 ) a_fit =  7.491671 ; b_fit =  2.466606 ;   R2 =  0.1648986 ;   RMSE =  8.95643 
## 8 ) a_fit =  14.61957 ; b_fit =  0.7389926 ;  R2 =  0.1518354 ;   RMSE =  8.95934 
## 8 ) a_fit =  36.13157 ; b_fit =  0.03005452 ;     R2 =  0.1783001 ;   RMSE =  9.1046 
## 8 ) a_fit =  109.8534 ; b_fit =  -0.01710957 ;    R2 =  0.2021991 ;   RMSE =  9.39526 
## 8 ) a_fit =  108.4526 ; b_fit =  -0.03030211 ;    R2 =  0.2090256 ;   RMSE =  9.358248 
## 9 ) a_fit =  8.986937 ; b_fit =  2.206359 ;   R2 =  0.1272976 ;   RMSE =  9.070061 
## 9 ) a_fit =  7.650902 ; b_fit =  2.440692 ;   R2 =  0.1586471 ;   RMSE =  9.021347 
## 9 ) a_fit =  14.87923 ; b_fit =  0.7285043 ;  R2 =  0.1469295 ;   RMSE =  9.019683 
## 9 ) a_fit =  37.4024 ; b_fit =  0.02766957 ;  R2 =  0.1659178 ;   RMSE =  9.36852 
## 9 ) a_fit =  99.99953 ; b_fit =  -0.01418411 ;    R2 =  0.1435705 ;   RMSE =  9.552597 
## 9 ) a_fit =  100.1835 ; b_fit =  -0.02569391 ;    R2 =  0.1505604 ;   RMSE =  9.516102 
## 10 ) a_fit =  11.57471 ; b_fit =  1.960843 ;  R2 =  0.1132833 ;   RMSE =  9.217313 
## 10 ) a_fit =  10.2941 ; b_fit =  2.127497 ;   R2 =  0.1355196 ;   RMSE =  9.063766 
## 10 ) a_fit =  17.07697 ; b_fit =  0.6752103 ;     R2 =  0.1347896 ;   RMSE =  9.052663 
## 10 ) a_fit =  43.64423 ; b_fit =  0.02357225 ;    R2 =  0.09105606 ;  RMSE =  9.576446 
## 10 ) a_fit =  80.67912 ; b_fit =  -0.006875877 ;  R2 =  0.05818177 ;  RMSE =  9.600197 
## 10 ) a_fit =  80.468 ; b_fit =  -0.01210274 ;     R2 =  0.05960173 ;  RMSE =  9.592897 
## 11 ) a_fit =  3.369919 ; b_fit =  3.344318 ;  R2 =  0.1408751 ;   RMSE =  9.132297 
## 11 ) a_fit =  3.217137 ; b_fit =  3.458154 ;  R2 =  0.2017893 ;   RMSE =  9.003591 
## 11 ) a_fit =  8.643776 ; b_fit =  1.013074 ;  R2 =  0.1740437 ;   RMSE =  9.007777 
## 11 ) a_fit =  41.28896 ; b_fit =  0.02641921 ;    R2 =  0.1010583 ;   RMSE =  9.485659 
## 11 ) a_fit =  105.4118 ; b_fit =  -0.01461405 ;   R2 =  0.1917603 ;   RMSE =  9.27535 
## 11 ) a_fit =  103.7901 ; b_fit =  -0.02526568 ;   R2 =  0.1948443 ;   RMSE =  9.245671 
## 12 ) a_fit =  2.8717 ; b_fit =  3.513259 ;    R2 =  0.153327 ;    RMSE =  9.169528 
## 12 ) a_fit =  4.033029 ; b_fit =  3.199581 ;  R2 =  0.1682395 ;   RMSE =  9.06045 
## 12 ) a_fit =  9.612162 ; b_fit =  0.9600938 ;     R2 =  0.154464 ;    RMSE =  9.040218 
## 12 ) a_fit =  40.85059 ; b_fit =  0.02620225 ;    R2 =  0.1183588 ;   RMSE =  9.714413 
## 12 ) a_fit =  92.04569 ; b_fit =  -0.01179612 ;   R2 =  0.1512497 ;   RMSE =  9.75852 
## 12 ) a_fit =  91.31656 ; b_fit =  -0.02042833 ;   R2 =  0.1523465 ;   RMSE =  9.732931 
## 13 ) a_fit =  13.97952 ; b_fit =  1.765124 ;  R2 =  0.113335 ;    RMSE =  9.107881 
## 13 ) a_fit =  12.26431 ; b_fit =  2.01955 ;   R2 =  0.1525834 ;   RMSE =  9.115751 
## 13 ) a_fit =  18.36556 ; b_fit =  0.6738635 ;     R2 =  0.1539626 ;   RMSE =  9.050074 
## 13 ) a_fit =  47.63929 ; b_fit =  0.01957556 ;    R2 =  0.08607894 ;  RMSE =  9.465863 
## 13 ) a_fit =  77.01419 ; b_fit =  -0.004729239 ;  R2 =  0.04437182 ;  RMSE =  9.686164 
## 13 ) a_fit =  77.0385 ; b_fit =  -0.008274721 ;   R2 =  0.04717384 ;  RMSE =  9.672191
colnames(r2df) <- c("Weeks", "a", "b", "R2", "RMSE", "VI")
num_cols <- c("Weeks", "a", "b", "R2", "RMSE")
r2df[num_cols] <- lapply(r2df[num_cols], as.numeric)
str(r2df)
## 'data.frame':    78 obs. of  6 variables:
##  $ Weeks: num  1 1 1 1 1 1 2 2 2 2 ...
##  $ a    : num  114.5 127.6 104 85.6 69.4 ...
##  $ b    : num  -1.7754 -1.8145 -0.8232 -0.1388 -0.0173 ...
##  $ R2   : num  0.1048 0.0724 0.0941 0.0534 0.085 ...
##  $ RMSE : num  9.37 9.42 9.35 9.45 9.97 ...
##  $ VI   : chr  "NDVI" "GNDVI" "EVI2" "SR" ...
r2df$Fieldname <- rep(selfield, length(lvl))
r2df$Approach <- rep("2b", length(lvl))
r2df_all <- rbind(r2df_all, r2df)

All vegetation indices + All weeks + All fields

load("Exponential_results_March31.RData", verbose = TRUE)
## Loading objects:
##   results_df

Plotting performance data - All VI + All weeks + All fields + All approaches

# "DMD_GH1", "SLS_ABH", "SLS_NS", "SSF_121" "SSF_66", "SSF_202", "PSF_111", "PSF_12"  
# selfield <- "DMD_GH1"

# results_df$yield_type <- recode_factor(factor(sel_dat4$yield_type), Pred_Yield = "Predicted yield", Yield_Mg.ha = "Actual yield")
  
results_df$VI <- factor(results_df$VI, levels = c("TGI", "EXG", "SR", "GNDVI", "EVI2", "NDVI"))

# selVI <- "EVI2"

sub_results <- results_df %>% 
  # filter(Fieldname == selfield) %>% 
  droplevels()

bar_plot<-ggplot(data=sub_results, aes(x=Weeks, y=R2, fill=VI)) +
  # geom_bar(stat="identity", position=position_dodge(0.8), width = 0.7) +
  geom_bar(stat="identity", position=position_dodge(), width = 0.6) +
  theme_bw() +
  theme(legend.position = "right", legend.box = "horizontal") +
  # theme(legend.position = "none") +
  # guides(colour = guide_legend(nrow = 1)) + 
  # ggtitle(selfield) + 
  # facet_wrap(~Approach) +
  facet_grid(Fieldname~Approach) + 
  # facet_grid(Fieldname~Approach) + 
  scale_x_continuous(breaks = seq(1, 14, 1)) +
  # scale_y_continuous(expand = c(0, 0)) + 
  scale_y_continuous(expand = c(0, 0), breaks = seq(0.0, 1.05, 0.2), limits = c(0.0, 1.05)) +
  # geom_text(aes(label=round(R2,2)), vjust=0.5, hjust = -0.2,  color="black",
  #           position = position_dodge(0), size=2, angle = 90)+
  # scale_fill_brewer(palette="Reds") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.title = element_text(size = 10, face = "bold"),
        # axis.title.y = element_blank(),
        axis.text = element_text(size = 8, color = "black"),
        # axis.text.y = element_blank(),
        # axis.ticks.length=unit(-1.5, "mm"),
        plot.title = element_text(size=14, face = "bold")) +
  theme(legend.title = element_text(face = "bold", size=8, margin = margin(0,0,7,0)),
        legend.text = element_text(vjust = 0.5),
        strip.text = element_text(face = "bold", size = 10)) 
  
bar_plot

# ggplotly(bar_plot)

# ggsave("All_VI_performance.pdf", width = 12, height = 8, units = "in")

Plotting performance data - NDVI

# "DMD_GH1", "SLS_ABH", "SLS_NS", "SSF_121" "SSF_66", "SSF_202", "PSF_111", "PSF_12"  
# selfield <- "SSF_66"

# results_df$yield_type <- recode_factor(factor(sel_dat4$yield_type), Pred_Yield = "Predicted yield", Yield_Mg.ha = "Actual yield")
  
results_df$Fieldname <- factor(results_df$Fieldname, levels = c("DMD_GH1", "SLS_ABH", "SLS_NS", "SSF_121", "SSF_66", "SSF_202", "PSF_111", "PSF_12"))

selVI <- "NDVI"

sub_results <- results_df %>% 
  filter(VI == selVI) %>% 
  droplevels()

bar_plot<-ggplot(data=sub_results, aes(x=Weeks, y=R2, fill=VI)) +
  geom_bar(stat="identity", position=position_dodge(0.8), width = 0.7, fill = "lightseagreen") +
  theme_bw() +
  # theme(legend.position = "bottom", legend.box = "horizontal") + 
  theme(legend.position = "none") +
  guides(colour = guide_legend(nrow = 1)) + 
  ggtitle(selVI) + 
  facet_grid(Fieldname~Approach) + 
  scale_x_continuous(breaks = seq(1, 14, 1)) +
  scale_y_continuous(expand = c(0, 0), breaks = seq(0.0, 1.05, 0.2), limits = c(0.0, 1.05)) +
  geom_text(aes(label=round(R2,2)), vjust=0.5, hjust = -0.2,  color="black",
            position = position_dodge(0), size=2, angle = 90)+
  # scale_fill_brewer(palette="Reds") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.title = element_text(size = 10, face = "bold"),
        # axis.title.y = element_blank(),
        axis.text = element_text(size = 8, color = "black"),
        # axis.text.y = element_blank(),
        # axis.ticks.length=unit(-1.5, "mm"),
        plot.title = element_text(size=14, face = "bold")) +
  theme(legend.title = element_text(face = "bold", size=8, margin = margin(0,0,7,0)),
        legend.text = element_text(vjust = 0.5),
        strip.text = element_text(face = "bold", size = 10)) 
  
bar_plot

# ggsave("VI_performance.pdf", width = 8.5, height = 11, units = "in")

Plotting performance data - EVI2

# "DMD_GH1", "SLS_ABH", "SLS_NS", "SSF_121" "SSF_66", "SSF_202", "PSF_111", "PSF_12"  
# selfield <- "SSF_66"

# results_df$yield_type <- recode_factor(factor(sel_dat4$yield_type), Pred_Yield = "Predicted yield", Yield_Mg.ha = "Actual yield")
  
results_df$Fieldname <- factor(results_df$Fieldname, levels = c("DMD_GH1", "SLS_ABH", "SLS_NS", "SSF_121", "SSF_66", "SSF_202", "PSF_111", "PSF_12"))

selVI <- "EVI2"

sub_results <- results_df %>% 
  filter(VI == selVI) %>% 
  droplevels()

bar_plot<-ggplot(data=sub_results, aes(x=Weeks, y=R2, fill=VI)) +
  geom_bar(stat="identity", position=position_dodge(0.8), width = 0.7, fill = "lightseagreen") +
  theme_bw() +
  # theme(legend.position = "bottom", legend.box = "horizontal") + 
  theme(legend.position = "none") +
  guides(colour = guide_legend(nrow = 1)) + 
  ggtitle(selVI) + 
  facet_grid(Fieldname~Approach) + 
  scale_x_continuous(breaks = seq(1, 14, 1)) +
  scale_y_continuous(expand = c(0, 0), breaks = seq(0.0, 1.05, 0.2), limits = c(0.0, 1.05)) +
  geom_text(aes(label=round(R2,2)), vjust=0.5, hjust = -0.2,  color="black",
            position = position_dodge(0), size=2, angle = 90)+
  # scale_fill_brewer(palette="Reds") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.title = element_text(size = 10, face = "bold"),
        # axis.title.y = element_blank(),
        axis.text = element_text(size = 8, color = "black"),
        # axis.text.y = element_blank(),
        # axis.ticks.length=unit(-1.5, "mm"),
        plot.title = element_text(size=14, face = "bold")) +
  theme(legend.title = element_text(face = "bold", size=8, margin = margin(0,0,7,0)),
        legend.text = element_text(vjust = 0.5),
        strip.text = element_text(face = "bold", size = 10)) 
  
bar_plot

# ggsave("VI_performance.pdf", width = 8.5, height = 11, units = "in")

Plotting performance data - GNDVI

# "DMD_GH1", "SLS_ABH", "SLS_NS", "SSF_121" "SSF_66", "SSF_202", "PSF_111", "PSF_12"  
# selfield <- "SSF_66"

# results_df$yield_type <- recode_factor(factor(sel_dat4$yield_type), Pred_Yield = "Predicted yield", Yield_Mg.ha = "Actual yield")
  
results_df$Fieldname <- factor(results_df$Fieldname, levels = c("DMD_GH1", "SLS_ABH", "SLS_NS", "SSF_121", "SSF_66", "SSF_202", "PSF_111", "PSF_12"))

selVI <- "GNDVI"

sub_results <- results_df %>% 
  filter(VI == selVI) %>% 
  droplevels()

bar_plot<-ggplot(data=sub_results, aes(x=Weeks, y=R2, fill=VI)) +
  geom_bar(stat="identity", position=position_dodge(0.8), width = 0.7, fill = "lightseagreen") +
  theme_bw() +
  # theme(legend.position = "bottom", legend.box = "horizontal") + 
  theme(legend.position = "none") +
  guides(colour = guide_legend(nrow = 1)) + 
  ggtitle(selVI) + 
  facet_grid(Fieldname~Approach) + 
  scale_x_continuous(breaks = seq(1, 14, 1)) +
  scale_y_continuous(expand = c(0, 0), breaks = seq(0.0, 1.05, 0.2), limits = c(0.0, 1.05)) +
  geom_text(aes(label=round(R2,2)), vjust=0.5, hjust = -0.2,  color="black",
            position = position_dodge(0), size=2, angle = 90)+
  # scale_fill_brewer(palette="Reds") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.title = element_text(size = 10, face = "bold"),
        # axis.title.y = element_blank(),
        axis.text = element_text(size = 8, color = "black"),
        # axis.text.y = element_blank(),
        # axis.ticks.length=unit(-1.5, "mm"),
        plot.title = element_text(size=14, face = "bold")) +
  theme(legend.title = element_text(face = "bold", size=8, margin = margin(0,0,7,0)),
        legend.text = element_text(vjust = 0.5),
        strip.text = element_text(face = "bold", size = 10)) 
  
bar_plot

# ggsave("VI_performance.pdf", width = 8.5, height = 11, units = "in")

Week 6 and 7 - Three VI + All approaches + All fields

sel_dat <- results_df %>% 
  filter(Weeks == 6 | Weeks == 7) %>% 
  filter(VI == "NDVI" | VI == "GNDVI" | VI == "EVI2") %>% 
  droplevels()

bar_plot<-ggplot(data=sel_dat, aes(x=Weeks, y=R2, fill=VI)) +
  geom_bar(stat="identity", position=position_dodge(0.8), width = 0.7) +
  theme_bw() +
  # theme(legend.position = "bottom", legend.box = "horizontal") + 
  theme(legend.position = "right") +
  guides(colour = guide_legend(nrow = 1)) + 
  # ggtitle(selVI) + 
  facet_grid(Fieldname~Approach) + 
  scale_x_continuous(breaks = seq(1, 14, 1)) +
  scale_y_continuous(expand = c(0, 0), breaks = seq(0.0, 1.05, 0.2), limits = c(0.0, 1.05)) +
  # geom_text(aes(label=round(R2,2)), vjust=0.5, hjust = -0.2,  color="black",
  #           position = position_dodge(0), size=2, angle = 90)+
  # scale_fill_brewer(palette="Reds") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.title = element_text(size = 10, face = "bold"),
        # axis.title.y = element_blank(),
        axis.text = element_text(size = 8, color = "black"),
        # axis.text.y = element_blank(),
        # axis.ticks.length=unit(-1.5, "mm"),
        plot.title = element_text(size=14, face = "bold")) +
  theme(legend.title = element_text(face = "bold", size=8, margin = margin(0,0,7,0)),
        legend.text = element_text(vjust = 0.5),
        strip.text = element_text(face = "bold", size = 10)) 
  
bar_plot

Saving the combined metrics data into CSV